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ABSTRACT 

Recent years have seen the discovery of many tidal streams through the Galaxy. Rela- 
tively straightforward observations of a stream allow one to deduce three phase-space 
coordinates of an orbit. An algorithm is presented that reconstructs the missing phase- 
space coordinates from these data. The reconstruction starts from assumed values of 
the Galactic potential and a distance to one point on the orbit, but with noise-free 
data the condition that energy be conserved on the orbit enables one to reject incorrect 
assumptions. The performance of the algorithm is investigated when errors are added 
to the input data that are comparable to those in published data for the streams of 
Pal 5. It is found that the algorithm returns distances and proper motions that are 
accurate to of order one percent, and enables one to reject quite reasonable but incor- 
rect trial potentials. In practical applications it will be important to minimize errors 
in the input data, and there is considerable scope for doing this. 
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1 INTRODUCTION 

Mapping the matter distribution in the outer reaches of 
the Galaxy, which is expected to be dark-matter domi- 
nated, is currently one of the central tasks in astrophysics. 
Recent surveys of the outer Galaxy have revealed nu- 
merous enhancements in the density of stars that have 
been secu rely identified as t i dal streamers from Galactic 
satellites ilOdenkirchen et al.l l200ll : iMaiewski et all |2004 
iBelokurov et al. 20061 ). Such structures effectively delineate 
the orbit in the Ga l actic force-field of the parent body 
(Johnston et al ] |l99fllOdenkirchen et af]l2003l ). Knowledge 
of the phase-space coordinates along an individual orbit pro- 
vides strong constraints on the force-field that confines the 
orbit, and thus on the matter distribution that gives rise 
to the force field. Consequently, in recent years consider- 
able effort has been devoted to the problem of determining 
the orbit of a satellite from observatio ns of its tidal stream 
|Law et alj|2005l : IFellhauer et~aT]|2007l ). 

Since tidal streams are generally discovered as enhance- 
ments in the density of stars on the sky, the phase-space 
coordinates of the satellite orbit that are easiest to deter- 
mine are its Galactic coordinates [l(u),b(u)], where u is an 
arbitrary parameter that increases along the orbit. With a 
moderate observational follow-up effort, it is always possible 
to augment these data with line-of-sight velocities v\\ (u). The 
other three coordinates (heliocentric distance s and proper 
motion components) are much harder to obtain, and are gen- 
erally either unknown or measured with very low accuracy. 



Given the incomplete nature of the data, the standard ap- 
proach to orbit fitting is to adopt a gravitational potential, 
and then to seek an orbit in this potentia l that is consistent 
with the data (e.g. IFellhauer et al.ll2007l ). The fitting pro- 
cedure is gene rally trial-and-error, a nd the final fit is often 
imperfect (e.g. iFellhauer et "al]|2006l ). 

Here we show that given the three coordinates (I, b, vu) 
along a section of an orbit in a gravitational potential $o(r), 
one can uniquely determine what the remaining three phase- 
space coordinates must be if the data are to be part of an 
orbit in a trial potential < I > t(r). If $t = $o, the recovered 
phase-space coordinates are consistent with conservation of 
energy, while when the trial and true potentials differ, en- 
ergy conservation appears to be violated. Therefore this pro- 
cedure can be used to hunt for the Galaxy's gravitational 
potential, and yields predictions for distances and proper 
motions along a tidal stream that can be tested observa- 
tionally. 



2 THE ALGORITHM 

We work in the inertial coordinate system in which the 
Galactic centre is at rest. Consequently, line-of-sight veloci- 
ties in this frame will differ from measured heliocentric veloc- 
ities by the projection of the Sun's velocity (whch we assume 
known) onto the line of sight. Let s be the position vector 
of the satellite with respect to an observer at the position of 
the Sun, and let F t (s) be the trial gravitational acceleration 
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(F t = — V$t)- Then if the trial potential is correct, we have 
s = F t , and vu — s.v — s so vu — F« + v.ds/dt, where Fm is 
the component of Ft along the line of sight. Also 

diss) .„ ds 
v = , — ss + s — 
dt dt 



(1) 



so 



(2) 



Moreover, the component of v in the plane of the sky vj 
satisfies 



|2 2-2 

\Vj_ I = s u 



s 2 (b 2 + cos 2 bf 



(3) 



where we have now fixed the meaning of the parameter u 
to be the angular distance along the trajectory. Eliminating 
|fj_| 2 between equations ((2} and ((3]), we obtain a quadratic 
equation for dt that has solution 



dt = 



2-Fii 



dv H ± A /dt) 2 - 4sFiidit 2 



(4) 



Thus if we guess s at some point on the trajectory, we can 
solve for dt between observations, and calculate E. We get 
a new value of E for each pair of data points. We can also 
update s from ds = vudt. 

To obtain a robust numerical scheme we divide ds = 
Wiidt and equation Q by the differential of u to obtain the 
set of coupled ordinary differential equations 



dt_ 

du 



1 

2F« 



dvu 
du 




isF« 



dt 



(5) 



ds 

d^ =W|l d^ 

The sign ambiguity in the first equation is resolved such that 
t increases along the trajectory. 



3 TESTS 

The algorithm was tested by reconstructing orbits from 
pseudodata obtained by projecting orbits in a model po- 
tential. The tests showed that it is important to determine 
accurately the distance u along the trajectory to each data 
point. This can be determined in a two-step process: we first 
obtain a crude estimate 



A * = E \l " b ^ 2 + cos2 [sfo+i + b Mh+i 



h) 2 - (6) 



Then we fit cubic splines in A to the values bi — b(\i) and 
li = i(Aj), and use the resulting functions of A to evaluate 
numerically the integrals 

r x i 

(7) 



dA yjb' 2 + cos 2 (b)l' 2 



where a prime denotes differentiation with respect to 
lambda. Finally we fit cubic splines to the values bi — b(ui), 
etc., and use these functions of u and the fourth-order 
Runge-Kutta algorithm to solve equations At each in- 
tegration step we re-evaluate E — \{v 2 \ + s 2 u 2 ) + <I>t. 

Orbits were integrated in the Miy amoto-Nagai poten- 
tial that has b/a = 0.2 and GM/a = 1 ( Miyamoto fc Nagail 
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Figure 1. Right panel: the projection onto the sky of a nu- 
merically integrated orbit in a Miyamoto-Nagai potential with 
b/a = 0.2. Left panel: the line-of-sight velocity as a function of b. 
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Figure 2. The log to base 10 of the rms variation in the energy 
when the orbit shown in Fig. [T] is reconstructed from an assumed 
initial distance s rather than its true value so. The triangles show 
results obtained when the reconstruction employs true potential, 
while the squares and pentagons are for less flattened potentials. 



1 19751 . and Fig.©. They were then reconstructed using either 
this potential or a Miyamoto-Nagai potential with different 
parameters. The 42 dots in the right panel of Fig. [T]show the 
input coordinates on the sky of the model orbit when it is 
viewed from a point in the symmetry plane and 8a from the 
centre, while the left panel shows the linbe-of-sight veloci- 
ties along the orbit. The orbit has energy — 0.0322GA//a. ft 
starts at b = —5 deg and a distance 15a from the point of 
observation, and moves out to a distance 22.6a as it rises to 
b = 68.9 deg. The triangles in Fig. [2] show the rms variation 
AE in the output energy E along the reconstructed orbit 
when the true potential is used to integrate equations ([5| 
from these data points and an initial distance s that is the 
stated multiple of its true value, so. There is a sharp min- 
imum in AE when s = so. The squares and pentagons in 
Fig. [2] show the corresponding results when the wrong po- 
tential is used to reconstruct the orbit - larger values of b/a 
generate rounder potentials. 

When the potential is correct, we find to a good ap- 
proximation that AE oc |1 — s/sq\, while with an incorrect 
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Figure 3. A contour plot of log 10 (Ai?) as a function of the pa- 
rameters used in the reconstruction of the orbit segment shown 
in Fig. [[] 




Figure 4. Isodensity contours in the meridional planes of 
Miyamoto-Nagai models with b/a = 0.2 (top) and b/a = 0.22 
(bottom). 



potential the curve of AE can be fitted near its minimum 
with AE oc - s/s ) 2 + k 2 with k / 0. Thus m a prac- 
tical case the quantity k can be used as a measure of how 
far the potential employed differs from its true value. 

Fig. [3]explores the effect of varying both the parameters 
b/a and GM/a of the potential used in the reconstruction of 
this orbit. Although an error in one parameter can to some 
extent be offset by an error in the other parameter, there is a 
sharp minimum in AE when the parameters take their true 
values. Thus from this one orbit one can determine both the 
flattening and the mass of the potential with precision. Fig. 2] 
illustrates how high this precision is by showing isodensity 
contours of the Miyamoto-Nagai models that have b/a = 0.2 
(top) and 0.22. Although these models are extremely similar, 
analysis of the orbit segment easily distinguishes them. 

Changing the length of the orbit segment analysed by 
a factor 3 either side of the length shown in Fig. [T] does 
not change the results significantly. In particular, the algo- 
rithm identifies incorrect potentials even with quite short 
orbit segments. 



These tests demonstrate that with high-quality data the 
algorithm has no difficulty recovering the true orbit if the 
potential is the correct one, and provides excellent discrim- 
ination against incorrect potentials. 



3.1 Multiple orbit segments 

Different streams will be sensitive to different aspects of the 
Galactic potential - streams near the plane will be sensitive 
to the structure of the disc, while the Magellanic Steam must 
be insensitive to the disc but sensitive to the total Galac- 
tic mass. The Galactic potential would be most strongly 
constrained by a plot like that shown in Fig. [3] in which 
the quantity contoured is the sum of AE for all available 
streams. 

Obscuration by dust near the plane sometimes breaks a 
stream into two or more segments with an unobserved gap 
between them. In principle each segment could be analysed 
independently, but this would not be the optimal procedure. 
First, all segments should be used to determine the path of 
the underlying orbit on the sky, including through obscured 
regions. Then the algorithm should be used to reconstruct 
each segment, and the rms energy variation along all recon- 
structed segments (treated as a whole) minimized within 
the space spanned by the initial distances assumed for each 
section. 



3.2 Effect of errors in the data 

In practice the input data will be less exact than in the tests 
above because tidal streams have finite widths on the sky 
and do not exactly trace a single orbit. In favourable cases 
they delineate an orbit with great precision: for example the 
tidal stream of a globular cluster such as that of Pal 5 is 
preci sely bounded by the cluster's orbit (jOdenkirchen et all 
120031 ). The tidal streams of larger bodies, such as the Sgr 
dwarf galaxy, are broad, and their relation to the orbit of 
the progenitor is less evident. Correspondingly, in the case 
of a globular cluster, the uncertainty in the radial velocity 
along the stream will be limited to measurement errors in the 
velocities of stars, while in the case of a dwarf galaxy one has 
to contend also with the velocity dispersion of the progen- 
itor, which both widens the range of velocities encountered 
at any point on the stream and leads to a systematic offset 
between the mean velocity in the stream and the velocity of 
the underlying orbit. 

When applying the algorithm to a stream such as that 
of Pal 5, one would not use the positions and velocities of 
individual stars directly but would make a judgment as to 
the relation of the orbit to the stellar stream, and estimate 
the radial velocity on the orbit by averaging the velocities of 
nearby stars and then making a correction for the expected 
offset between the stream and orbital velocities. When all 
this had been done, one would contrive that the trajectory 
of the orbit on the sky and in radial velocity was smooth. 
Therefore, in tests it is not appropriate to ad independent 
errors to each data point, but rather to add correlated errors 
that leave the trajectory through (/, b, vu) space smooth. We 
do this as follows. For x = v«,b,l cos 6, we add a Fourier 
series 
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Figure 5. As for Fig.[2]except that equation JS} has been used to 
add errors to b and I cos(fe) that each amount to 2 arcmin rms. The 
minimum value of AE for b/a = 0.5 is bigger than for b/a = 0.2 
by a factor 1.3. 
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Figure 7. Data points for an orbit segment ten times as long as 
in Fig. [1] 
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Figure 6. As for Fig. [2] except that equation (JSj has been used 
to add errors to Dm that amount to 0.1 percent of the circular 
speed at the start of the segment. The minimum value of AE for 
b/a = 0.5 is bigger than for b/a = 0.2 by a factor 1.1. 



Sx(u) — N X n cos(2ri7m + f n ), 



(8) 



where n = 0—9 is an integer, u = u/u max ranges from zero to 
unity along the trajectory, and £ n is randomly distributed in 
(0,2-7r). The amplitudes X n determine the power spectrum 
of the noise. We have no way of knowing what this would 
be, but have assumed that X n oc (1 + n) -11 . 

Studies of the tidal stream of Pal 5 provide a guide to 
the magnitude of the errors to be encou ntered in real data. 
The stream's full width is ~ 30 arcmin |Odenkirchen et al.l 
2003), so it should be possible to determine the position 
of the underlying orbit with an error that does not ex- 
ceed a tenth of this, an d could be substantially smaller. 
lOdenkirchen et ail (|2002l ) report velocities along the stream 
that have errors ~ 0.15 kms" 1 plus a dispersion of order 
0.9 km s -1 arising from binarity. The latter dispersion could 
be beaten down by making observations at several epochs; 
even at one epoch, measurements of many stars would en- 
able the velocity of the underlying orbit to be obtained to 
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Figure 8. The results of reconstructing 100 realizations of the 
orbit segment of Fig. [7] with astrometric errors of 2 arcmin and 
velocity errors of 0.1 percent of v c . The points show the values of 
s/sq and b/a that gave the smallest value of AE, together with 
that value of AE. 



a precision of at least 0.2 km s^ 1 . Hence we focus here on 
errors of 2 arcmin in b and Zcos(&), and errors in vu that are 
0.1 percent of the circular speed at the start of the orbit 
segment. 

Fig. [5] shows that with these errors in b and I cos(fe) 
it is still possible to determine the distance to the orbit 
with exquisite accuracy, and to tell that the potential with 
b/a = 0.2 is to be preferred to that with b/a = 0.5. 

Fig.[6]shows the results obtained when the radial veloc- 
ities are in error by 0.1 percent of the circular speed at the 
starting point of the segment. The impact of these velocity 
errors is larger than that of the position errors. However, it 
is still just possible to reject the potential with b/a — 0.5, 
and the distance is accurately determined even if the wrong 
potential is chosen. When the orbit segment used for the 
reconstructions covers a time interval that is longer by a 
factor ten (Fig. [7]) , it becomes possible to reject the poten- 
tial with b/a — 0.5 with security. Fig. [8] illustrates this find- 
ing by showing for each of 100 different realizations of the 
noise added to the positions and velocities of Fig.[7]the best 
fitting values of s/so and b/a together with the associated 
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Figure 9. Variation in reconstructed energy E and time T along 
the orbit shown in Fig. [7] when the correct potential and distance 
are used when v\\ is in error by 0.1 percent of v c . The vertical 
dotted lines mark the times when the orbit passes through the 
plane. 



values of AE. All distance estimates lie within 2 percent of 
the truth; the mean of s/sq = 1.002 while its rms dispersion 
is 0.012. Values of b/a in the range (0.1, 1) were searched 
but no value bigger that 0.35 was recovered; the mean of 
b/a = 0.201 while its rms dispersion was 0.095. Thus with 
these errors in the data it is possible to say that the potential 
is highly flattened, but not actually determe b/a. 

Fig. [5] shows the reconstructed values of time and en- 
ergy along the best reconstruction of the long segment plot- 
ted in Fig. [7] when the velocities are in error by 0.1 pecent 
of the circular speed. Vertical dotted lines mark the instants 
at which the orbit passes through the plane; pericentre and 
the point of closest approach to the Sun lie between pairs 
of these plane-crossings. Some combination of plane cross- 
ing and the peak in the proper motion as the orbit passes 
closest to the Sun induces features in both the time and en- 
ergy plots. These features are independent of the accuracy 
parameter specified for the integration of equations ([5]), and 
they have no analogues when error-free data are used for the 
reconstruction. The reconstructed time T is systematically 
later than it should be by a few percent, and the recon- 
structed energy is correspondingly larger than it should be. 

Since the proper motions that the algorithm predicts are 
inversely proportional to the elapse of T along the stream, 
the figure shows that even with errors in the input velocities, 
proper motions are generally accurate to a few percent. 

The algorithm's sensitivity to errors in v\\ arises because 
the underlying physical idea is that the change in vu between 
data points requires a certain time interval, and this time 
interval determines both the change in distance between the 
points and the orbit's proper motion. Thus the reconstruc- 
tion is driven by differences in Vu , which are sensitive to 
errors in vu . When equation ([8| is used to simulate errors, 
the error in the difference in v\\ between neighbouring points 
scales like u^ixi so longer orbit segments produce more ac- 
curate results. 



4 CONCLUSIONS 

A tidal stream delineates the orbit of its progenitor, so it can 
be used to obtain three phase-space coordinates of an orbit. 
To complete the specification of the orbit, one needs to asso- 
ciate with each of its points on the sky a distance and a time 
of passage. The algorithm presented here yields these addi- 
tional quantities for an assumed potential and a distance to 
one point on the stream. Unless these two inputs are correct, 
there is no guarantee that energy will be conserved along the 
reconstructed orbit, and indeed numerical experiments show 
that with noise-free data, the rms variation in energy along 
the reconstructed orbit allows one to identify the potential 
and determine the correct initial distance. 

The algorithm is less sensitive to realistic errors in the 
position of the stream on the sky than to likely errors in 
U|l because the reconstruction is driven by differences in vu 
between grid points. 

The algorithm predicts the distances and proper mo- 
tions of objects in the stream, so when these data are avail- 
able, they will provide an interesting check on the correct- 
ness of all inputs and assumptions. Crude distance estimates 
are generally available and seriously erroneous potentials can 
probably be excluded because they return astrophysically 
unacceptable distances. With presently attainable accuracy 
of the input data it should be possible to determine distances 
to with about one percent, and to constrain the Galactic po- 
tential quite strongly. 

It is important to recognize that the inputs to the al- 
gorithm should be the sky coordinates and line-of-sight ve- 
locities of the orbit of a test particle, not the correspond- 
ing coordinates of a stream; streams delineate orbits, but 
they do not actually lie along them. In practice a significant 
effort will be required to extract the optimum input data 
from measurements of any given stream, and an optimal ap- 
proach to this task does not appear to have been published. 
However, the experiments presented here indicate that the 
rewards for obtaining accurate input data will be very high: 
even one orbit segment provides a wealth of precision in- 
formation about the Galactic potential and observationally 
testable predictions about the stream. The Galactic poten- 
tial could be even more strongly constrained by analyzing 
several orbit segments together. 
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